Load Model Data

Description of two models:

Two simple models are implemented in ACT-R, performing the gambling task depending on two distinctive mechanisms.

Model1: Declarative Model (version 2.0)

The model randomly evaluates an action GUESS-MORE or GUESS-LESS, then retrieving memories about prior feedback. If “WIN” memory is retrieved, the model chooses the evaluated action. If “LOSE”/“Neutral” memory is retrieved, the model chooses the alternative action.

Model2: RL Models.

The model makes decision relying on procedural memory, pressing “MORE”/“LESS” key. GUESS-MORE and GUESS-LESS productions are competing. In the end, a feedback is given, the model receives either +1/-1/0 reward to all previous productions, and the utility of all previous fired productions will be affected.

model1.raw <- read.csv("./MODEL120210424_100307_fnca_log.csv") #%>% rename(Trial=X)
model2.raw <- read.csv("./MODEL220210424_100307_fnca_log.csv") #%>% rename(Trial=X)

#model1 <- read.csv("MODEL120210420.csv")
#model2 <- read.csv("MODEL220210420.csv")

model1.raw <- model1.raw %>% 
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2),
         RT = as.numeric(RT),
         Epoch = as.numeric(Epoch)) 

model2.raw <- model2.raw %>%
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2),
         RT = as.numeric(RT),
         Epoch = as.numeric(Epoch)) 

# only 
#model1 <- model1.raw %>% filter(ans==.1 & bll==.1 & lf==.1)
#model2 <- model2.raw %>% filter(alpha==.1 & egs==.1 & r==.1)

Compute win-stay probabilities

future_moves <- function(responses) {
  c(responses[2:length(responses)], NA)
}
past_moves <- function(responses) {
  c(NA, responses[1:length(responses)-1])
}

count_responses <- function(model.dat) {
  model.dat <- model.dat %>%
   mutate(FutureResponse = future_moves(CurrentResponse),
         PastResponse = past_moves(CurrentResponse),
         PreviousFeedback = past_moves(TrialType),
         ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1))
  return(model.dat)
}

# NA for all first trial in each block
clean_previous_future <- function(model.dat) {
  model.dat <- model.dat %>%
    mutate(PastResponse=ifelse(BlockTrial==0, NA, PastResponse),
           PreviousFeedback=ifelse(BlockTrial==0, NA, PreviousFeedback),
           FutureResponse=ifelse(BlockTrial==7, NA, FutureResponse),
           ResponseSwitch=ifelse(BlockTrial==7, NA, ResponseSwitch),
           )
  return(model.dat)
}

#model1.count <- clean_previous_future(count_responses(model1.raw))

Aggregate PSwitch, RT by ParamID, TrialType, BlockType

### This aggregate looks at overall, what is the average pattern of 50 simulation (50 epoch) across different parameter set

# aggregate data by parameter set, BlockType, TrialType
m1.aggregate.ByParam.CurrentResponse.BlockType <- function(m1.dat) {
  res <- m1.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
    group_by(ans, bll, lf, BlockType, TrialType) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT)) %>%
    ungroup() %>%
    mutate(ParamID = group_indices(., ans, bll, lf))
  return(res)
}

m2.aggregate.ByParam.CurrentResponse.BlockType <- function(m2.dat) {
  res <- m2.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
    group_by(egs, alpha, r, BlockType, TrialType) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT)) %>%
    ungroup() %>%
    mutate(ParamID = group_indices(., egs, alpha, r))
  return(res)
}

# aggregate data by parameter set, BlockType, PreviousFeedback
m1.aggregate.ByParam.PreviousFeedback.BlockType <- function(m1.dat) {
  res <- m1.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
    group_by(ans, bll, lf, BlockType, PreviousFeedback) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))%>%
    ungroup() %>%
    mutate(ParamID = group_indices(., ans, bll, lf))
  return(res)
}

m2.aggregate.ByParam.PreviousFeedback.BlockType <- function(m2.dat) {
  res <- m2.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
    group_by(egs, alpha, r, BlockType, PreviousFeedback) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))%>%
    ungroup() %>%
    mutate(ParamID = group_indices(., egs, alpha, r))
  return(res)
}


#model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)

Aggregate PSwitch, RT by Epoch, TrialType, BlockType. (This only fit for one param set)

### This aggregate looks at for single parameter set, how consistent 50 runs are

# aggregate data by Epoch, BlockType, PreviousFeedback
aggregate.ByEpoch.CurrentResponse.BlockType <- function(m.dat) {
  res <- m.dat %>%
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
  group_by(Epoch, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
  return (res)
}

aggregate.ByEpoch.PreviousFeedback.BlockType <- function(m.dat) {
  res <- m.dat %>%
  filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
  group_by(Epoch, BlockType, PreviousFeedback) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
  return (res)
}
#model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)
# plot subject aggregated data, PSwitch by TrialType x BlockType
subj.PSwitch.aggplot <- function(subj.dat, subjID) {
  res <- ggplot(subj.dat %>% filter(HCPID==subjID), 
                  aes(x = TrialType, y = PSwitch)) +
    facet_grid( ~ BlockType) +
    geom_line(aes(group = BlockType), alpha=.5) + 
    geom_point(alpha=.5, size=3, aes(col=TrialType)) +
    ylim(0,1) +
    theme_pander()
  return (res)
}

# plot model aggregated data, PSwitch by TrialType x BlockType x ParamID/Epoch
m.PSwitch.aggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = TrialType, y = PSwitch, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(m.dat, aes(x = TrialType, y = PSwitch, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } else {
    print("Choose one: ParamID or Epoch")
  }
  return (res)
}

# plot model aggregated data, T by TrialType x BlockType x ParamID/Epoch
m.RT.aggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = TrialType, y = RT, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(m.dat, aes(x = TrialType, y = RT, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } else {
    print("Choose one: ParamID or Epoch")
  }
  return (res)
}


#m.PSwitch.aggplot(model1.aggregate, "ParamID")
#m.RT.aggplot(model1.aggregate, "ParamID")

Let’s look at the PSwitch and RT as a function of Current TrialType

model1.count <- clean_previous_future(count_responses(model1.raw))
model2.count <- clean_previous_future(count_responses(model2.raw))


model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)
model2.aggregate <- m2.aggregate.ByParam.CurrentResponse.BlockType(model2.count)
m.PSwitch.aggplot(model1.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

m.PSwitch.aggplot(model2.aggregate, "ParamID") + ggtitle("MODEL2")
## [1] "Choose one: ParamID or Epoch"

m.RT.aggplot(model1.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

m.RT.aggplot(model2.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

Summary of Pswitch vs. CurrentTrial

Model1:

Model2:

Summary of RT vs. CurrentTrial

Model1:

Model2:

model1.pf_aggregate <- m1.aggregate.ByParam.PreviousFeedback.BlockType(model1.count)
model2.pf_aggregate <- m2.aggregate.ByParam.PreviousFeedback.BlockType(model2.count)
# plot subject aggregated data, PSwitch by PreviousFeedback x BlockType
subj.PSwitch.pfaggplot <- function(subj.dat, subjID) {
  res <- ggplot(subj.dat %>% filter(HCPID==subjID), 
                  aes(x = PreviousFeedback, y = PSwitch)) +
    facet_grid( ~ BlockType) +
    geom_line(aes(group = BlockType), alpha=.5) + 
    geom_point(alpha=.5, size=3, aes(col=PreviousFeedback)) +
    ylim(0,1) +
    theme_pander()
  return (res)
}

# plot model aggregated data, PSwitch by PreviousFeedback x BlockType x ParamID/Epoch
m.PSwitch.pfaggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = PreviousFeedback, y = PSwitch, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(model1.aggregate, aes(x = PreviousFeedback, y = PSwitch, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } else {
    print("Choose one: ParamID or Epoch")
  }
  return (res)
}

# plot model aggregated data, T by TrialType x BlockType x ParamID/Epoch
m.RT.pfaggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = PreviousFeedback, y = RT, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(model1.aggregate, aes(x = PreviousFeedback, y = RT, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } else {
    print("Choose one: ParamID or Epoch")
  }
  return (res)
}

Let’s visualize it

m.PSwitch.pfaggplot(model1.pf_aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

m.PSwitch.pfaggplot(model2.pf_aggregate, "ParamID") + ggtitle("MODEL2")
## [1] "Choose one: ParamID or Epoch"

aov(PSwitch ~ BlockType * PreviousFeedback, model2.pf_aggregate  %>% filter(PreviousFeedback!="Neutral"))  %>%
  anova_summary() %>%
  xtable() %>%
  kable(digits=4) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

aov(RT ~ BlockType * PreviousFeedback, model2.pf_aggregate  %>% filter(PreviousFeedback!="Neutral"))  %>%
  anova_summary() %>%
  xtable() %>%
  kable(digits=4) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Summary of Pswicth vs. PreviousTrial

Model1:

Model2:

Summary of RT vs. PreviousTrial

Model1:

Model2:

blockId <- sort(rep(1:4, 8))

gambling <- gambling %>% 
  group_by(HCPID, RunNumber) %>%
  mutate(BlockId = blockId, 
         Response = if_else(is.na(QuestionMark.RESP), 0, QuestionMark.RESP))

bg <- gambling %>%
  group_by(BlockId, HCPID, RunNumber) %>%
  summarise(BlockType = paste("Mostly", Mode(TrialType)))

gambling <- inner_join(gambling, bg)

past_moves <- function(responses) {
  c(NA, responses[1:length(responses)-1])
}

gambling <- gambling %>%
  group_by(RunNumber, BlockId) %>%
  mutate(CurrentResponse = Response,
         FutureResponse = future_moves(Response),
         PastResponse = past_moves(Response),
         PreviousFeedback = past_moves(TrialType),
         RT = QuestionMark.RT/1200) %>%
  filter(FutureResponse != 0) %>%
  #filter(CurrentResponse != 0) %>%
  mutate(ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1))   %>%
  select(HCPID, RunNumber, BlockId, RT, Response, PastResponse,
         CurrentResponse, FutureResponse, ResponseSwitch, BlockType, TrialType,
         PreviousFeedback, FeedbackImage)

gambling$PreviousFeedback <- as_factor(gambling$PreviousFeedback)
gambling$PreviousFeedback <- revalue(gambling$PreviousFeedback,
                                         c("1"="Reward", 
                                           "2"="Punishment", 
                                           "3"="Neutral"))
subj1.aggregate <- aggregate.CurrentTrial(gambling)
subj1.pf_aggregate <- aggregate.PreviousTrial(gambling)

PSwitch.aggplot(subj1.aggregate) + ggtitle("SUBJ - 100307")
RT.aggplot(subj1.aggregate) + ggtitle("SUBJ - 100307")

PSwitch.pfaggplot(subj1.pf_aggregate) + ggtitle("SUBJ - 100307")
RT.pfaggplot(subj1.pf_aggregate) + ggtitle("SUBJ - 100307")
gambling <- read_tsv("../../gambling_data.txt")
gambling.trials <- gambling %>% 
  #filter(RunNumber==1) %>%
  select(HCPID, Trial, RunNumber, TrialType, RunTrialNumber, Block, 
         QuestionMark.RESP, QuestionMark.ACC, QuestionMark.RT,
         R1MostlyReward1, R1MostlyReward2, R1MostlyPunishment1, R1MostlyPunishment3, QuestionMark.OnsetDelay, 
         ConsecSameResp, ConsecLargerGuesses, ConsecSmallerGuesses, ConsecRTLess200) %>% 
  rename(CurrentResponse=QuestionMark.RESP, RT=QuestionMark.RT)  %>%
  mutate(FutureResponse = future_moves(CurrentResponse),
         PastResponse = past_moves(CurrentResponse), 
         ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1),
         PreviousFeedback = past_moves(as.character(TrialType)),
         BlockType = case_when(
           !is.na(R1MostlyReward1) ~ "MostlyReward", 
           !is.na(R1MostlyReward2) ~ "MostlyReward", 
           !is.na(R1MostlyPunishment1) ~ "MostlyPunishment", 
           !is.na(R1MostlyPunishment3) ~ "MostlyPunishment"
         )) %>%
  select(-R1MostlyReward1, -R1MostlyReward2, -R1MostlyPunishment1, -R1MostlyPunishment3) %>%
  filter(!is.na(TrialType))

#write.csv(gambling.trials, "../../gambling_clean_data.csv")
gambling.aggregate <- gambling.trials.session1 %>% 
  filter(!is.na(ResponseSwitch) & !is.na(RT)) %>%
  group_by(HCPID, TrialType, BlockType, PreviousFeedback) %>%
  summarise(PSwitch = mean(ResponseSwitch),RT = mean(RT))

ggplot(gambling.aggregate,
       aes(x = TrialType, y = PSwitch, col = TrialType)) +
  facet_grid( ~ BlockType, labeller = label_both) +
  geom_point(position = position_jitter(width = 0.1, height = 0.05),
             alpha = 0.1) +
  scale_color_brewer(palette = "Set2") +
  stat_summary(fun.data = "mean_cl_boot", col="black") +
  theme_pander()

gambling.trials <- gambling.trials.session1 %>% select(HCPID, Trial, RunTrialNumber, TrialType, BlockType)

gambling.trials.session1 %>%
  gghistogram("RT", add = "median")

Load Subject Data

gambling.clean <- read_csv("../../bin/gambling_clean_data.csv") %>% 
  select(-BlockType) %>% rename(BlockType=BlockTypeCoded)


gambling.cfagg <- gambling.clean %>%
  na.omit() %>%
  group_by(HCPID, BlockType, TrialType) %>%
  summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

gambling.pfagg <- gambling.clean %>%
  na.omit() %>%
  group_by(HCPID, BlockType, PreviousFeedback) %>%
  summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

# long format
gambling.cfagg <- left_join(gambling.cfagg %>% select(-RT) %>%
  spread(key = TrialType, value = PSwitch, sep = "."),
  gambling.cfagg %>% select(-PSwitch) %>%
  spread(key = TrialType, value = RT, sep = "."),
  by = c("HCPID", "BlockType"), suffix=c(".PSwitch",".RT"))

gambling.pfagg <- left_join(gambling.pfagg %>% select(-RT) %>%
  spread(key = PreviousFeedback, value = PSwitch, sep = "."),
  gambling.pfagg %>% select(-PSwitch) %>%
  spread(key = PreviousFeedback, value = RT, sep = "."),
  by = c("HCPID", "BlockType"), suffix=c(".PSwitch",".RT"))

#write_csv(gambling.cfagg, "../../bin/gambling_cfagg.csv")
#write_csv(gambling.pfagg, "../../bin/gambling_pfagg.csv")

Load clean data and select one subj data

Note: here, ResponseSwitch is 1 if: CurrentResponse != FutureResponse; is 0 if CurrentResponse==FutureResponse

gambling.clean <- read_csv("../../bin/gambling_clean_data.csv") %>% mutate(BlockType=BlockTypeCoded)
subjID = "100307_fnca"

#  Pswitch as a function of current trial type and block type
subj.cfaggregate = gambling.clean %>%
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>% 
  group_by(HCPID, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch, rm.na=T), RT = mean(RT, rm.na=T))

#  Pswitch as a function of previous trial type and block type
subj.pfaggregate = gambling.clean %>%
  filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback)) %>% 
  group_by(HCPID, BlockType, PreviousFeedback) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch, rm.na=T), RT = mean(RT, rm.na=T))
# test whether prop is calculate correctly
subj.cfaggregate %>% filter(HCPID==subjID)
gambling.clean %>% filter(HCPID==subjID) %>%
  filter(!is.na(ResponseSwitch)) %>%
  group_by(BlockType, TrialType , ResponseSwitch) %>%
  dplyr::summarise(n=n()) %>%
  mutate(freq = n / sum(n)) %>%
  xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))

#>> yes correct

Visualize subject Pswitch pattern as a function of previous feedback

subj.PSwitch.aggplot(subj.cfaggregate, subjID) + ggtitle(paste("SUB", subjID, sep = "-"))

subj.PSwitch.pfaggplot(subj.pfaggregate, subjID) + ggtitle(paste("SUB", subjID, sep = "-"))

Calculate RMSE and Find Best Fit Parameter

Distribution of model1 parameters :ans, :bll, :lf

Note: This is not grid-search parameter fitting. The parameter space was chosen by optimization function by adding null_penalty (count of null responses)

dim(model1.raw)
## [1] 150400     11
model1.full.count <- clean_previous_future(count_responses(model1.raw))

model1.full.agg <- model1.full.count %>% 
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
  group_by(ans, bll, lf, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

model1.full.agg %>% gghistogram(x="ans", binwidth = .05)

model1.full.agg %>% gghistogram(x="bll", binwidth = .05)

model1.full.agg %>% gghistogram(x="lf", binwidth = .05)

Distribution model2 parameter :egs, :alpha, and :r

dim(model2.raw)
## [1] 339200     11
model2.full.count <- clean_previous_future(count_responses(model2.raw))

model2.full.agg <- model2.full.count %>% 
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
  group_by(egs, alpha, r, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

model2.full.agg %>% gghistogram(x="egs", binwidth = .05)

model2.full.agg %>% gghistogram(x="alpha", binwidth = .05)

model2.full.agg %>% gghistogram(x="r", binwidth = .05)

Next, we calalculate RMSE for each parameter combination, find the best parameter with min RMSE

find_best_parameter <- function(model, model.agg, sub.agg) {
  if (model=="model1") {
    m1.param <- unique(model.agg[,c('ans','bll','lf')])
    #m1.param=m1.param[1:10,]
    m1.param$RMSE = NA
    m1.param$RMSE = as.numeric(m1.param$RMSE)
    
    for (i in 1:nrow(m1.param)) {
      m1.agg <- model.agg %>% 
        filter(ans==m1.param[[i,'ans']] & bll==m1.param[[i,'bll']] & lf==m1.param[[i,'lf']])
        rmse_value <- rmse(m1.agg$PSwitch, sub.agg$PSwitch)
      m1.param[[i, 'RMSE']] <- rmse_value
    }
    return (m1.param)
  } else {
    m2.param <- unique(model.agg[,c('egs','alpha','r')])
    #m2.param=m2.param[1:10,]
    m2.param$RMSE = NA
    m2.param$RMSE = as.numeric(m2.param$RMSE)
    
    for (i in 1:nrow(m2.param)) {
      m2.agg <- model.agg %>% 
        filter(egs==m2.param[[i,'egs']] & alpha==m2.param[[i,'alpha']] & r==m2.param[[i,'r']])
        rmse_value <- rmse(m2.agg$PSwitch, sub.agg$PSwitch)
      m2.param[[i, 'RMSE']] <- rmse_value
    }
    return (m2.param)
  }
}

m1.bestp <- find_best_parameter("model1", model1.full.agg, subj.cfaggregate)
m2.bestp <- find_best_parameter("model2", model2.full.agg, subj.cfaggregate)
m1.bestp$model = "model1"
m2.bestp$model = "model2"
m1.bestp <- m1.bestp[order(m1.bestp$RMSE),]
m2.bestp <- m2.bestp[order(m2.bestp$RMSE),]

m1.bestp %>% head() %>% xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))
ans bll lf RMSE model
1.2265 0.3263 0.757 0.2738 model1
1.2257 0.1000 0.100 0.2744 model1
1.2380 0.1000 0.100 0.2753 model1
1.2265 0.3263 0.739 0.2754 model1
3.0902 0.1000 0.100 0.2758 model1
1.2265 0.3263 0.729 0.2770 model1 
m2.bestp %>% head() %>% xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))
egs alpha r RMSE model
0.7295 0.1038 0.6972 0.2730 model2
0.7295 0.1038 2.5030 0.2730 model2
2.1291 0.1122 5.4124 0.2732 model2
2.1287 0.1122 7.5503 0.2738 model2
0.7295 0.1038 1.4486 0.2739 model2
0.7295 0.1044 0.1000 0.2740 model2 
m1.bestp[[1,'ans']]
## [1] 1.226471

By looking at min RMSE:

Best parameter chosen for model1 :ans=1.2264712 :bll=0.3262827, :lf=0.757

Best parameter chosen for model2 is :egs=0.7294902 ; alpha=0.1037721, :r=0.6972462

Load optimal simulation

Having the optimal parameter set, we can run simulation 800 times again by setting optimal parameter to the model.

subj1 = "./no_seed/MODEL220210425_100307_fnca_best800.csv"
subj2 = "./seed/MODEL220210425_100307_fnca_best800.csv"
subjID = "100307_fnca"
model2.best800 <- read.csv(subj1) %>%
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2),
         RT = as.numeric(RT),
         Epoch = as.numeric(Epoch)) 
model2.best800.count <- clean_previous_future(count_responses(model2.best800))
model2.best800.aggregate <- aggregate.ByEpoch.CurrentResponse.BlockType(model2.best800.count)


ggarrange(m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==0), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==1), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==2), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==3), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==4), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==5), "Epoch"))

m.PSwitch.aggplot(model2.best800.aggregate, "Epoch")

Calculate Log-Likelihood

model2.best800.aggregate <- model2.best800.aggregate %>%
  group_by(BlockType, TrialType) %>%
  dplyr::summarise(PSwitch.m = mean(PSwitch), PSwitch.sd = sd(PSwitch))

subj1.aggregate <- subj.cfaggregate %>% filter(HCPID == subjID) %>% rename(PSwitch.subj = PSwitch)

# bic 
bic.ll <- function(k, n, logL) {
  return(k * log(n) - 2 * logL)
}
#= n log(RSS/n) + (p + 1) log(n)
bic.rmse <- function(k, n, rss) {
  return(n * log(rss/n) - (k + 1) * log(n))
}

model2.best800.LL <- left_join(model2.best800.aggregate, subj1.aggregate) %>%
  mutate(PSwitch.z = (PSwitch.m-PSwitch.subj)/(PSwitch.sd),
         PSwitch.probz = dnorm(PSwitch.z), 
         PSwitch.logprobz = log(PSwitch.probz)) %>%
  group_by(HCPID) %>%
  dplyr::summarise(PSwitch.LL = sum(PSwitch.logprobz), RMSE = rmse(PSwitch.m, PSwitch.subj), RSS = RMSE^2 * 6) 

model2.best800.LL %>% mutate(BIC.LL = bic.ll(k=3, n=6, PSwitch.LL), BIC.rmse = bic.rmse(k=3, n=6, RSS)) %>%
  xtable() %>%
  kable(digits=4) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
HCPID PSwitch.LL RMSE RSS BIC.LL BIC.rmse
100307_fnca -8.1648 0.1824 0.1996 21.7049 -27.5868